#!/usr/bin/env julia

"""
    test_crtm_interface.jl

Comprehensive test suite for CRTM (Community Radiative Transfer Model) interface in GSI.
Tests radiative transfer computations, Jacobian calculations, surface emissivity,
and integration with GSI observation operators.

Test Coverage:
- CRTM initialization and configuration
- Forward radiative transfer calculations
- Jacobian (K-matrix) computations for atmospheric profiles
- Surface emissivity modeling for different surface types
- Channel-specific radiance simulations
- Integration with satellite observation operators
- Performance testing for operational usage
"""

using Test
using GSICoreAnalysis
using GSICoreAnalysis.ObservationProcessing.CRTMInterface
using LinearAlgebra
using Statistics
using Random
using Dates

"""Test CRTM initialization and configuration"""
@testset "CRTM Initialization" begin
    println("  Testing CRTM initialization...")
    
    @testset "Sensor Configuration" begin
        # Test AMSU-A sensor configuration
        amsu_config = CRTMSensorConfig(
            sensor_id = "amsu_a",
            satellite_id = "noaa_18",
            n_channels = 15,
            channel_list = collect(1:15)
        )
        
        @test amsu_config.sensor_id == "amsu_a"
        @test amsu_config.satellite_id == "noaa_18"
        @test amsu_config.n_channels == 15
        @test length(amsu_config.channel_list) == 15
        
        # Test MHS sensor configuration
        mhs_config = CRTMSensorConfig(
            sensor_id = "mhs",
            satellite_id = "noaa_19",
            n_channels = 5,
            channel_list = [1, 2, 3, 4, 5]
        )
        
        @test mhs_config.sensor_id == "mhs"
        @test mhs_config.n_channels == 5
        
        # Test IASI sensor configuration (partial channels)
        iasi_channels = [16, 38, 49, 51, 55, 57, 59, 61, 63, 66, 
                        70, 72, 74, 79, 81, 83, 85, 87, 89, 92]
        iasi_config = CRTMSensorConfig(
            sensor_id = "iasi",
            satellite_id = "metop_a",
            n_channels = length(iasi_channels),
            channel_list = iasi_channels
        )
        
        @test iasi_config.n_channels == 20
        @test iasi_config.channel_list == iasi_channels
    end
    
    @testset "Atmospheric Profile Setup" begin
        # Create standard atmosphere profile
        n_levels = 42  # Typical GSI vertical levels
        
        profile = CRTMAtmosphericProfile{Float64}(
            n_levels = n_levels,
            pressure = zeros(n_levels),
            temperature = zeros(n_levels),
            water_vapor = zeros(n_levels),
            ozone = zeros(n_levels)
        )
        
        # US Standard Atmosphere
        pressure_levels = reverse([
            1013.25, 898.76, 795.01, 701.21, 616.60, 540.48, 471.18,
            407.00, 348.25, 294.30, 244.58, 199.53, 158.65, 121.39,
            87.20, 55.51, 26.71, 11.97, 5.75, 2.87, 1.49, 0.80,
            0.43, 0.23, 0.12, 0.065, 0.035, 0.018, 0.010, 0.005,
            0.003, 0.002, 0.001, 0.0007, 0.0005, 0.0003, 0.0002,
            0.0001, 0.00007, 0.00005, 0.00003, 0.00001
        ])
        
        profile.pressure .= pressure_levels .* 100.0  # Convert to Pa
        
        # Temperature profile
        for k in 1:n_levels
            p = pressure_levels[k]
            if p > 226.32  # Troposphere
                profile.temperature[k] = 288.15 - 6.5 * (1013.25 - p) / 100.0
            elseif p > 54.75  # Stratosphere
                profile.temperature[k] = 216.65
            else  # Mesosphere
                profile.temperature[k] = 216.65 + 2.8 * (54.75 - p) / 10.0
            end
        end
        
        # Water vapor profile (exponential decay)
        for k in 1:n_levels
            p = pressure_levels[k]
            profile.water_vapor[k] = 0.01 * exp(-0.0001 * (1013.25 - p))
        end
        
        # Ozone profile (peak in stratosphere)
        for k in 1:n_levels
            p = pressure_levels[k]
            if 10.0 <= p <= 100.0  # Stratospheric ozone peak
                profile.ozone[k] = 8e-6 * exp(-((log(p) - log(30.0))/0.8)^2)
            else
                profile.ozone[k] = 1e-7
            end
        end
        
        @test profile.n_levels == n_levels
        @test length(profile.pressure) == n_levels
        @test length(profile.temperature) == n_levels
        @test all(profile.temperature .> 180.0)  # Reasonable temperatures
        @test all(profile.temperature .< 320.0)
        @test all(profile.water_vapor .>= 0.0)   # Non-negative humidity
        @test all(profile.ozone .>= 0.0)         # Non-negative ozone
    end
    
    @testset "Surface Property Setup" begin
        # Test various surface types
        
        # Ocean surface
        ocean_surface = CRTMSurfaceProperties{Float64}(
            surface_type = SURFACE_OCEAN,
            land_fraction = 0.0,
            water_fraction = 1.0,
            ice_fraction = 0.0,
            snow_fraction = 0.0,
            surface_temperature = 285.0,
            wind_speed = 7.5,
            wind_direction = 225.0,
            salinity = 35.0
        )
        
        @test ocean_surface.surface_type == SURFACE_OCEAN
        @test ocean_surface.water_fraction == 1.0
        @test ocean_surface.salinity == 35.0
        
        # Land surface
        land_surface = CRTMSurfaceProperties{Float64}(
            surface_type = SURFACE_LAND,
            land_fraction = 1.0,
            water_fraction = 0.0,
            ice_fraction = 0.0,
            snow_fraction = 0.0,
            surface_temperature = 295.0,
            vegetation_fraction = 0.8,
            soil_moisture = 0.25,
            leaf_area_index = 3.5
        )
        
        @test land_surface.surface_type == SURFACE_LAND
        @test land_surface.land_fraction == 1.0
        @test land_surface.vegetation_fraction == 0.8
        
        # Sea ice surface
        ice_surface = CRTMSurfaceProperties{Float64}(
            surface_type = SURFACE_ICE,
            land_fraction = 0.0,
            water_fraction = 0.0,
            ice_fraction = 1.0,
            snow_fraction = 0.3,
            surface_temperature = 265.0,
            ice_thickness = 1.5,
            snow_depth = 0.15
        )
        
        @test ice_surface.surface_type == SURFACE_ICE
        @test ice_surface.ice_fraction == 1.0
        @test ice_surface.surface_temperature < 273.15  # Below freezing
    end
end

"""Test forward radiative transfer calculations"""
@testset "Forward Radiative Transfer" begin
    println("  Testing forward radiative transfer...")
    
    @testset "AMSU-A Radiance Calculation" begin
        # Set up AMSU-A sensor
        amsu_config = CRTMSensorConfig(
            sensor_id = "amsu_a",
            satellite_id = "noaa_18",
            n_channels = 15,
            channel_list = collect(1:15)
        )
        
        # Create atmospheric profile
        n_levels = 42
        profile = create_standard_atmosphere(n_levels)
        
        # Ocean surface
        surface = CRTMSurfaceProperties{Float64}(
            surface_type = SURFACE_OCEAN,
            water_fraction = 1.0,
            surface_temperature = 285.0,
            wind_speed = 5.0
        )
        
        # Viewing geometry
        geometry = CRTMGeometry{Float64}(
            satellite_zenith_angle = 25.0,
            satellite_azimuth_angle = 180.0,
            solar_zenith_angle = 45.0,
            solar_azimuth_angle = 0.0
        )
        
        # Compute radiances
        radiances = crtm_forward(amsu_config, profile, surface, geometry)
        
        @test length(radiances) == 15
        @test all(radiances .> 180.0)   # Reasonable brightness temperatures
        @test all(radiances .< 320.0)
        
        # Test channel-specific properties
        # Channel 1 (23.8 GHz) - surface sensitive
        @test radiances[1] > 250.0      # Should be warm (surface influenced)
        
        # Channel 9 (57.29 GHz) - peak near 800 hPa
        @test 220.0 < radiances[9] < 280.0
        
        # Channel 14 (57.29 GHz) - peak near 300 hPa  
        @test 200.0 < radiances[14] < 250.0  # Cooler upper atmosphere
        
        # Channels should generally decrease with increasing peak altitude
        @test radiances[1] > radiances[5] > radiances[14]
    end
    
    @testset "MHS Water Vapor Channels" begin
        # Set up MHS sensor
        mhs_config = CRTMSensorConfig(
            sensor_id = "mhs",
            satellite_id = "noaa_19",
            n_channels = 5,
            channel_list = [1, 2, 3, 4, 5]
        )
        
        # Create moist atmospheric profile
        profile = create_standard_atmosphere(42)
        
        # Increase water vapor content
        profile.water_vapor .*= 3.0  # Triple the humidity
        
        ocean_surface = CRTMSurfaceProperties{Float64}(
            surface_type = SURFACE_OCEAN,
            water_fraction = 1.0,
            surface_temperature = 288.0,
            wind_speed = 8.0
        )
        
        geometry = CRTMGeometry{Float64}(
            satellite_zenith_angle = 20.0,
            satellite_azimuth_angle = 90.0,
            solar_zenith_angle = 30.0
        )
        
        # Compute radiances for moist profile
        radiances_moist = crtm_forward(mhs_config, profile, ocean_surface, geometry)
        
        @test length(radiances_moist) == 5
        
        # Create dry atmospheric profile for comparison
        profile_dry = deepcopy(profile)
        profile_dry.water_vapor .*= 0.3  # Reduce humidity to 30%
        
        radiances_dry = crtm_forward(mhs_config, profile_dry, ocean_surface, geometry)
        
        # Water vapor channels should show sensitivity to humidity changes
        # Channel 3 (183.31±1 GHz) - most sensitive to water vapor
        @test abs(radiances_moist[3] - radiances_dry[3]) > 5.0  # Significant difference
        
        # Channel 1 (89.0 GHz) - window channel, less sensitive
        @test abs(radiances_moist[1] - radiances_dry[1]) < 10.0  # Smaller difference
        
        # Moist atmosphere generally produces lower brightness temperatures
        # in water vapor channels due to increased absorption
        @test radiances_moist[3] < radiances_dry[3]
    end
    
    @testset "Surface Type Sensitivity" begin
        # Common sensor and profile setup
        amsu_config = CRTMSensorConfig(
            sensor_id = "amsu_a",
            satellite_id = "noaa_18",
            n_channels = 5,
            channel_list = [1, 2, 6, 9, 14]  # Sample channels
        )
        
        profile = create_standard_atmosphere(42)
        
        geometry = CRTMGeometry{Float64}(
            satellite_zenith_angle = 15.0,
            satellite_azimuth_angle = 0.0
        )
        
        # Ocean surface
        ocean = CRTMSurfaceProperties{Float64}(
            surface_type = SURFACE_OCEAN,
            water_fraction = 1.0,
            surface_temperature = 285.0,
            wind_speed = 6.0
        )
        
        # Land surface  
        land = CRTMSurfaceProperties{Float64}(
            surface_type = SURFACE_LAND,
            land_fraction = 1.0,
            surface_temperature = 295.0,
            vegetation_fraction = 0.7,
            soil_moisture = 0.3
        )
        
        # Sea ice surface
        ice = CRTMSurfaceProperties{Float64}(
            surface_type = SURFACE_ICE,
            ice_fraction = 1.0,
            surface_temperature = 260.0,
            ice_thickness = 2.0
        )
        
        # Compute radiances for each surface type
        radiances_ocean = crtm_forward(amsu_config, profile, ocean, geometry)
        radiances_land = crtm_forward(amsu_config, profile, land, geometry)
        radiances_ice = crtm_forward(amsu_config, profile, ice, geometry)
        
        # Surface-sensitive channels should show significant differences
        # Channel 1 (23.8 GHz) - highly surface sensitive
        ocean_land_diff = abs(radiances_ocean[1] - radiances_land[1])
        ocean_ice_diff = abs(radiances_ocean[1] - radiances_ice[1])
        
        @test ocean_land_diff > 5.0   # Significant surface sensitivity
        @test ocean_ice_diff > 10.0   # Ice should be quite different
        
        # Upper atmospheric channels should be less sensitive
        # Channel 5 (higher peaking) - less surface sensitive
        upper_ocean_land_diff = abs(radiances_ocean[5] - radiances_land[5])
        @test upper_ocean_land_diff < ocean_land_diff  # Less surface sensitivity
        
        # Ice surface should generally be colder
        @test radiances_ice[1] < radiances_ocean[1]
        @test radiances_ice[1] < radiances_land[1]
    end
end

"""Test Jacobian calculations"""
@testset "Jacobian Calculations" begin
    println("  Testing Jacobian calculations...")
    
    @testset "Temperature Jacobians" begin
        # Setup for Jacobian testing
        amsu_config = CRTMSensorConfig(
            sensor_id = "amsu_a",
            satellite_id = "noaa_18",
            n_channels = 10,
            channel_list = collect(1:10)
        )
        
        n_levels = 42
        profile = create_standard_atmosphere(n_levels)
        
        surface = CRTMSurfaceProperties{Float64}(
            surface_type = SURFACE_OCEAN,
            water_fraction = 1.0,
            surface_temperature = 285.0,
            wind_speed = 5.0
        )
        
        geometry = CRTMGeometry{Float64}(
            satellite_zenith_angle = 20.0,
            satellite_azimuth_angle = 0.0
        )
        
        # Compute temperature Jacobians
        temp_jacobians = crtm_jacobian_temperature(amsu_config, profile, surface, geometry)
        
        @test size(temp_jacobians) == (10, n_levels)  # n_channels × n_levels
        
        # Test Jacobian properties
        for ch in 1:10
            # Jacobians should sum to approximately 1 (unit response)
            jacobian_sum = sum(temp_jacobians[ch, :])
            @test 0.8 <= jacobian_sum <= 1.2
            
            # Peak sensitivity location should make sense
            peak_level = argmax(abs.(temp_jacobians[ch, :]))
            
            if ch <= 3  # Surface channels
                @test peak_level <= 10  # Peak near surface
            elseif ch >= 8  # Upper atmosphere channels
                @test peak_level >= 25  # Peak in upper levels
            end
        end
        
        # Test finite difference accuracy
        δT = 0.1  # Small temperature perturbation
        level_test = 20  # Mid-troposphere level
        
        # Reference radiances
        radiances_ref = crtm_forward(amsu_config, profile, surface, geometry)
        
        # Perturb temperature at test level
        profile_pert = deepcopy(profile)
        profile_pert.temperature[level_test] += δT
        
        radiances_pert = crtm_forward(amsu_config, profile_pert, surface, geometry)
        
        # Finite difference Jacobian
        fd_jacobians = (radiances_pert .- radiances_ref) ./ δT
        
        # Compare with analytical Jacobians
        for ch in 1:10
            analytical_jac = temp_jacobians[ch, level_test]
            fd_jac = fd_jacobians[ch]
            
            relative_error = abs(analytical_jac - fd_jac) / abs(fd_jac)
            @test relative_error < 0.05  # 5% accuracy
        end
    end
    
    @testset "Water Vapor Jacobians" begin
        # Setup for water vapor sensitivity
        mhs_config = CRTMSensorConfig(
            sensor_id = "mhs",
            satellite_id = "noaa_19", 
            n_channels = 5,
            channel_list = [1, 2, 3, 4, 5]
        )
        
        n_levels = 42
        profile = create_standard_atmosphere(n_levels)
        
        surface = CRTMSurfaceProperties{Float64}(
            surface_type = SURFACE_OCEAN,
            water_fraction = 1.0,
            surface_temperature = 288.0
        )
        
        geometry = CRTMGeometry{Float64}(
            satellite_zenith_angle = 25.0
        )
        
        # Compute water vapor Jacobians
        wv_jacobians = crtm_jacobian_water_vapor(mhs_config, profile, surface, geometry)
        
        @test size(wv_jacobians) == (5, n_levels)
        
        # Water vapor channels should have significant sensitivity
        # Channel 3 (183.31±1 GHz) - strong WV absorption
        ch3_max_jac = maximum(abs.(wv_jacobians[3, :]))
        ch1_max_jac = maximum(abs.(wv_jacobians[1, :]))  # Window channel
        
        @test ch3_max_jac > 2.0 * ch1_max_jac  # WV channel more sensitive
        
        # Peak sensitivity should be in troposphere for most channels
        for ch in [3, 4, 5]  # Water vapor channels
            peak_level = argmax(abs.(wv_jacobians[ch, :]))
            @test peak_level <= 30  # Tropospheric peak
        end
        
        # Test finite difference validation
        δq = 0.0001  # Small humidity perturbation (kg/kg)
        test_level = 15
        
        radiances_ref = crtm_forward(mhs_config, profile, surface, geometry)
        
        profile_pert = deepcopy(profile)
        profile_pert.water_vapor[test_level] += δq
        
        radiances_pert = crtm_forward(mhs_config, profile_pert, surface, geometry)
        
        fd_wv_jac = (radiances_pert .- radiances_ref) ./ δq
        
        for ch in 1:5
            analytical_jac = wv_jacobians[ch, test_level]
            fd_jac = fd_wv_jac[ch]
            
            if abs(fd_jac) > 1e-6  # Only test significant Jacobians
                relative_error = abs(analytical_jac - fd_jac) / abs(fd_jac)
                @test relative_error < 0.1  # 10% accuracy (WV Jacobians noisier)
            end
        end
    end
    
    @testset "Surface Parameter Jacobians" begin
        # Test surface temperature and emissivity Jacobians
        amsu_config = CRTMSensorConfig(
            sensor_id = "amsu_a",
            satellite_id = "noaa_18",
            n_channels = 6,
            channel_list = [1, 2, 3, 4, 15, 6]  # Mix of surface-sensitive channels
        )
        
        profile = create_standard_atmosphere(42)
        
        surface = CRTMSurfaceProperties{Float64}(
            surface_type = SURFACE_LAND,
            land_fraction = 1.0,
            surface_temperature = 295.0,
            vegetation_fraction = 0.6,
            soil_moisture = 0.25
        )
        
        geometry = CRTMGeometry{Float64}(
            satellite_zenith_angle = 30.0
        )
        
        # Compute surface temperature Jacobians
        tskin_jacobians = crtm_jacobian_surface_temperature(amsu_config, profile, surface, geometry)
        
        @test length(tskin_jacobians) == 6
        
        # Surface channels should have larger sensitivity
        # Channel 1 (23.8 GHz) - very surface sensitive
        @test tskin_jacobians[1] > 0.5  # Strong positive sensitivity
        
        # Higher peaking channels should have lower sensitivity
        if length(tskin_jacobians) >= 6
            @test tskin_jacobians[6] < tskin_jacobians[1]  # Less surface sensitive
        end
        
        # Test surface emissivity Jacobians (land surface)
        emissivity_jacobians = crtm_jacobian_surface_emissivity(amsu_config, profile, surface, geometry)
        
        @test length(emissivity_jacobians) == 6
        
        # All emissivity Jacobians should be positive for surface channels
        for ch in 1:3  # Surface-sensitive channels
            @test emissivity_jacobians[ch] > 0.0
        end
        
        # Validate with finite differences
        δTs = 0.5  # 0.5 K surface temperature perturbation
        
        radiances_ref = crtm_forward(amsu_config, profile, surface, geometry)
        
        surface_pert = deepcopy(surface)
        surface_pert.surface_temperature += δTs
        
        radiances_pert = crtm_forward(amsu_config, profile, surface_pert, geometry)
        
        fd_tskin_jac = (radiances_pert .- radiances_ref) ./ δTs
        
        for ch in 1:6
            analytical_jac = tskin_jacobians[ch]
            fd_jac = fd_tskin_jac[ch]
            
            relative_error = abs(analytical_jac - fd_jac) / abs(fd_jac)
            @test relative_error < 0.05  # 5% accuracy
        end
    end
end

"""Test surface emissivity modeling"""
@testset "Surface Emissivity Modeling" begin
    println("  Testing surface emissivity modeling...")
    
    @testset "Ocean Surface Emissivity" begin
        # Test ocean surface emissivity calculation
        frequencies = [23.8, 31.4, 50.3, 52.8, 89.0]  # GHz (AMSU-A channels)
        
        # Calm ocean conditions
        calm_ocean = CRTMSurfaceProperties{Float64}(
            surface_type = SURFACE_OCEAN,
            water_fraction = 1.0,
            surface_temperature = 285.0,
            wind_speed = 2.0,  # Calm
            salinity = 35.0
        )
        
        geometry = CRTMGeometry{Float64}(
            satellite_zenith_angle = 20.0
        )
        
        emissivity_calm = compute_ocean_emissivity(frequencies, calm_ocean, geometry)
        
        @test length(emissivity_calm) == length(frequencies)
        @test all(0.3 .<= emissivity_calm .<= 1.0)  # Physical range
        
        # Rough ocean conditions
        rough_ocean = CRTMSurfaceProperties{Float64}(
            surface_type = SURFACE_OCEAN,
            water_fraction = 1.0,
            surface_temperature = 285.0,
            wind_speed = 15.0,  # Rough
            salinity = 35.0
        )
        
        emissivity_rough = compute_ocean_emissivity(frequencies, rough_ocean, geometry)
        
        # Rough surfaces should have higher emissivity
        for i in 1:length(frequencies)
            @test emissivity_rough[i] >= emissivity_calm[i]
        end
        
        # Lower frequencies should generally have lower emissivity
        @test emissivity_calm[1] < emissivity_calm[end]  # 23.8 GHz < 89 GHz
        
        # Test angle dependence
        geometry_nadir = CRTMGeometry{Float64}(satellite_zenith_angle = 0.0)
        geometry_slant = CRTMGeometry{Float64}(satellite_zenith_angle = 50.0)
        
        emissivity_nadir = compute_ocean_emissivity(frequencies, calm_ocean, geometry_nadir)
        emissivity_slant = compute_ocean_emissivity(frequencies, calm_ocean, geometry_slant)
        
        # Emissivity should generally increase with viewing angle
        for i in 1:length(frequencies)
            @test emissivity_slant[i] >= emissivity_nadir[i] - 0.1  # Allow small tolerance
        end
    end
    
    @testset "Land Surface Emissivity" begin
        # Test land surface emissivity for different surface types
        frequencies = [23.8, 31.4, 50.3, 89.0, 150.0]
        
        # Desert surface
        desert = CRTMSurfaceProperties{Float64}(
            surface_type = SURFACE_LAND,
            land_fraction = 1.0,
            surface_temperature = 305.0,
            vegetation_fraction = 0.05,  # Sparse vegetation
            soil_moisture = 0.05,        # Very dry
            sand_fraction = 0.9
        )
        
        geometry = CRTMGeometry{Float64}(satellite_zenith_angle = 25.0)
        
        emissivity_desert = compute_land_emissivity(frequencies, desert, geometry)
        
        @test length(emissivity_desert) == length(frequencies)
        @test all(0.7 .<= emissivity_desert .<= 1.0)  # Dry land typically high emissivity
        
        # Forest surface
        forest = CRTMSurfaceProperties{Float64}(
            surface_type = SURFACE_LAND,
            land_fraction = 1.0,
            surface_temperature = 285.0,
            vegetation_fraction = 0.95,  # Dense vegetation
            soil_moisture = 0.35,        # Moist soil
            leaf_area_index = 5.0
        )
        
        emissivity_forest = compute_land_emissivity(frequencies, forest, geometry)
        
        # Forest should have different emissivity characteristics
        @test all(0.8 .<= emissivity_forest .<= 1.0)
        
        # Vegetation effects should be frequency dependent
        # Higher frequencies more sensitive to vegetation
        veg_effect_low = emissivity_forest[1] - emissivity_desert[1]   # 23.8 GHz
        veg_effect_high = emissivity_forest[end] - emissivity_desert[end]  # 150 GHz
        
        # This relationship can vary, but there should be differences
        @test abs(veg_effect_low - veg_effect_high) > 0.02
        
        # Wetland surface (high soil moisture)
        wetland = CRTMSurfaceProperties{Float64}(
            surface_type = SURFACE_LAND,
            land_fraction = 1.0,
            surface_temperature = 280.0,
            vegetation_fraction = 0.7,
            soil_moisture = 0.6,      # Very moist
            water_fraction = 0.3      # Some standing water
        )
        
        emissivity_wetland = compute_land_emissivity(frequencies, wetland, geometry)
        
        # Wet surfaces should have different characteristics
        @test all(0.6 .<= emissivity_wetland .<= 1.0)
        
        # Lower frequencies should show more water influence
        @test emissivity_wetland[1] < emissivity_desert[1]  # Water reduces emissivity at low freq
    end
    
    @testset "Snow and Ice Emissivity" begin
        frequencies = [23.8, 31.4, 50.3, 89.0, 150.0]
        
        # Fresh snow surface
        fresh_snow = CRTMSurfaceProperties{Float64}(
            surface_type = SURFACE_LAND,
            land_fraction = 1.0,
            surface_temperature = 268.0,
            snow_fraction = 1.0,
            snow_depth = 0.3,         # 30 cm fresh snow
            snow_density = 100.0,     # kg/m³ (fresh snow)
            snow_grain_size = 0.1     # mm (fine grains)
        )
        
        geometry = CRTMGeometry{Float64}(satellite_zenith_angle = 20.0)
        
        emissivity_fresh_snow = compute_snow_emissivity(frequencies, fresh_snow, geometry)
        
        @test length(emissivity_fresh_snow) == length(frequencies)
        @test all(0.5 .<= emissivity_fresh_snow .<= 1.0)
        
        # Fresh snow should have high emissivity at low frequencies
        @test emissivity_fresh_snow[1] > 0.9  # 23.8 GHz
        
        # Emissivity should decrease with frequency for snow
        @test emissivity_fresh_snow[1] > emissivity_fresh_snow[end]
        
        # Aged/wet snow
        aged_snow = CRTMSurfaceProperties{Float64}(
            surface_type = SURFACE_LAND,
            land_fraction = 1.0,
            surface_temperature = 271.0,  # Near melting
            snow_fraction = 1.0,
            snow_depth = 0.15,            # Compacted
            snow_density = 300.0,         # Denser
            snow_grain_size = 1.0,        # Larger grains
            liquid_water_content = 0.05   # Some liquid water
        )
        
        emissivity_aged_snow = compute_snow_emissivity(frequencies, aged_snow, geometry)
        
        # Aged snow should have lower emissivity, especially at high frequencies
        @test emissivity_aged_snow[end] < emissivity_fresh_snow[end]
        
        # Sea ice surface
        sea_ice = CRTMSurfaceProperties{Float64}(
            surface_type = SURFACE_ICE,
            ice_fraction = 1.0,
            surface_temperature = 263.0,
            ice_thickness = 2.0,
            snow_fraction = 0.2,      # Some snow cover
            snow_depth = 0.05
        )
        
        emissivity_ice = compute_ice_emissivity(frequencies, sea_ice, geometry)
        
        @test all(0.4 .<= emissivity_ice .<= 1.0)
        
        # Sea ice emissivity characteristics
        # Should be different from snow due to ice properties
        @test abs(mean(emissivity_ice) - mean(emissivity_fresh_snow)) > 0.05
    end
end

"""Test integration with GSI observation operators"""
@testset "GSI Integration" begin
    println("  Testing GSI integration...")
    
    @testset "Radiance Observation Operator Integration" begin
        # Create GSI configuration
        config = AnalysisConfig{Float64}(
            grid_size = (6, 6, 42),  # Match CRTM levels
            ensemble_size = 0
        )
        
        # Initialize state vector with realistic atmospheric state
        sv = StateVector(config)
        fill!(sv.values, 0.0)
        
        # Extract variable slices
        temp_slice = get_variable_slice(sv, "temperature")
        q_slice = get_variable_slice(sv, "specific_humidity")
        ps_slice = get_variable_slice(sv, "surface_pressure")
        
        # Initialize with standard atmosphere
        fill!(ps_slice, 101325.0)
        
        nx, ny, nz = config.grid_size
        for k in 1:nz
            for j in 1:ny
                for i in 1:nx
                    idx = linear_index(i, j, k, config.grid_size)
                    
                    # Height-based profiles
                    p_level = 101325.0 * exp(-0.2 * (k-1))  # Approximate pressure
                    
                    if p_level > 22632.0
                        temp_slice[idx] = 288.15 * (p_level/101325.0)^0.1903
                    else
                        temp_slice[idx] = 216.65
                    end
                    
                    q_slice[idx] = 0.015 * exp(-0.0001 * (101325.0 - p_level))
                end
            end
        end
        
        # Create satellite observation
        sat_obs = SatelliteObservation{Float64}(
            obs_type = "amsu_a",
            channel = 6,
            value = 242.5,
            error = 1.2,
            latitude = 40.0,
            longitude = -100.0,
            satellite_zenith_angle = 20.0,
            satellite_azimuth_angle = 180.0,
            scan_position = 15,
            time = DateTime(2024, 1, 1, 12, 0, 0)
        )
        
        # Apply integrated observation operator (should use CRTM)
        h_x = apply_satellite_radiance_operator(sat_obs, sv, config)
        
        @test 200.0 <= h_x <= 280.0  # Reasonable brightness temperature
        @test isfinite(h_x)
        
        # Test that different channels give different results
        sat_obs_ch2 = SatelliteObservation{Float64}(
            obs_type = "amsu_a",
            channel = 2,
            value = 255.0,
            error = 1.5,
            latitude = 40.0,
            longitude = -100.0,
            satellite_zenith_angle = 20.0,
            time = DateTime(2024, 1, 1, 12, 0, 0)
        )
        
        h_x_ch2 = apply_satellite_radiance_operator(sat_obs_ch2, sv, config)
        @test abs(h_x - h_x_ch2) > 2.0  # Different channels should differ significantly
    end
    
    @testset "Tangent Linear Integration" begin
        config = AnalysisConfig{Float64}(
            grid_size = (4, 4, 20),
            ensemble_size = 0
        )
        
        sv = StateVector(config)
        cv = ControlVector(config)
        
        # Initialize vectors
        randn!(sv.values)
        randn!(cv.values)
        
        # Scale to reasonable atmospheric values
        temp_slice = get_variable_slice(sv, "temperature")
        temp_slice .= temp_slice .* 20.0 .+ 250.0  # 230-270 K range
        
        q_slice = get_variable_slice(sv, "specific_humidity")
        q_slice .= abs.(q_slice) .* 0.005  # 0-0.01 kg/kg range
        
        sat_obs = SatelliteObservation{Float64}(
            obs_type = "amsu_a",
            channel = 8,
            value = 235.0,
            error = 1.8,
            latitude = 40.0,
            longitude = -100.0,
            satellite_zenith_angle = 25.0,
            time = DateTime(2024, 1, 1, 12, 0, 0)
        )
        
        # Test tangent linear operator with CRTM integration
        h_tl = apply_tangent_linear_radiance_operator(sat_obs, cv, sv, config)
        
        @test isfinite(h_tl)
        @test abs(h_tl) < 50.0  # Reasonable TL response
        
        # Test adjoint operator
        cv_adj = ControlVector(config)
        fill!(cv_adj.values, 0.0)
        
        innovation = 3.5  # obs - H(x)
        apply_adjoint_radiance_operator!(sat_obs, innovation, cv_adj, sv, config)
        
        # Adjoint should modify control vector
        @test norm_cv(cv_adj) > 0.0
        
        # Test adjoint property (approximately)
        left_side = h_tl * innovation
        right_side = dot_product(cv, cv_adj)
        
        relative_error = abs(left_side - right_side) / abs(right_side)
        @test relative_error < 0.05  # 5% accuracy for adjoint property
    end
end

"""Test performance characteristics"""
@testset "CRTM Performance" begin
    println("  Testing CRTM performance...")
    
    @testset "Computational Performance" begin
        # Test computational efficiency for operational use
        
        # Set up realistic configuration
        amsu_config = CRTMSensorConfig(
            sensor_id = "amsu_a",
            satellite_id = "noaa_18",
            n_channels = 15,
            channel_list = collect(1:15)
        )
        
        profile = create_standard_atmosphere(42)
        
        surface = CRTMSurfaceProperties{Float64}(
            surface_type = SURFACE_OCEAN,
            water_fraction = 1.0,
            surface_temperature = 285.0,
            wind_speed = 7.0
        )
        
        geometry = CRTMGeometry{Float64}(
            satellite_zenith_angle = 20.0,
            satellite_azimuth_angle = 180.0
        )
        
        # Time forward calculations
        n_calls = 100
        println("    Timing $(n_calls) CRTM forward calls...")
        
        forward_time = @elapsed begin
            for i in 1:n_calls
                radiances = crtm_forward(amsu_config, profile, surface, geometry)
            end
        end
        
        @test forward_time < 5.0  # Should complete within 5 seconds
        
        avg_time_per_call = forward_time / n_calls
        @test avg_time_per_call < 0.05  # Less than 50ms per call
        
        println("    Forward calculation time: $(forward_time:.2f) seconds total")
        println("    Average per call: $(avg_time_per_call*1000:.1f) ms")
        
        # Time Jacobian calculations
        n_jacobian_calls = 20  # Fewer due to higher computational cost
        println("    Timing $(n_jacobian_calls) CRTM Jacobian calls...")
        
        jacobian_time = @elapsed begin
            for i in 1:n_jacobian_calls
                temp_jac = crtm_jacobian_temperature(amsu_config, profile, surface, geometry)
                wv_jac = crtm_jacobian_water_vapor(amsu_config, profile, surface, geometry)
            end
        end
        
        @test jacobian_time < 10.0  # Should complete within 10 seconds
        
        avg_jacobian_time = jacobian_time / n_jacobian_calls
        @test avg_jacobian_time < 0.5  # Less than 500ms per call
        
        println("    Jacobian calculation time: $(jacobian_time:.2f) seconds total")
        println("    Average per call: $(avg_jacobian_time*1000:.1f) ms")
    end
    
    @testset "Memory Usage" begin
        # Test memory efficiency
        initial_memory = Base.gc_live_bytes()
        
        # Process multiple profiles
        n_profiles = 50
        amsu_config = CRTMSensorConfig(
            sensor_id = "amsu_a",
            satellite_id = "noaa_18",
            n_channels = 15,
            channel_list = collect(1:15)
        )
        
        for i in 1:n_profiles
            # Create varied profiles
            profile = create_standard_atmosphere(42)
            profile.temperature .+= randn(42) * 5.0  # Add variability
            profile.water_vapor .*= (1.0 + 0.5 * rand())  # Scale humidity
            
            surface = CRTMSurfaceProperties{Float64}(
                surface_type = rand() > 0.5 ? SURFACE_OCEAN : SURFACE_LAND,
                surface_temperature = 280.0 + 20.0 * rand(),
                wind_speed = 10.0 * rand()
            )
            
            geometry = CRTMGeometry{Float64}(
                satellite_zenith_angle = 60.0 * rand()
            )
            
            # Forward calculation
            radiances = crtm_forward(amsu_config, profile, surface, geometry)
            
            # Jacobian calculation (every 10th profile)
            if i % 10 == 0
                temp_jac = crtm_jacobian_temperature(amsu_config, profile, surface, geometry)
            end
            
            # Periodic garbage collection
            if i % 20 == 0
                GC.gc()
            end
        end
        
        final_memory = Base.gc_live_bytes()
        memory_growth = final_memory - initial_memory
        
        # Memory growth should be reasonable
        @test memory_growth < 100 * 1024 * 1024  # Less than 100 MB growth
        
        println("    Memory growth: $(memory_growth ÷ (1024*1024)) MB")
    end
end

"""Helper functions"""

function create_standard_atmosphere(n_levels::Int)
    """Create a standard atmosphere profile for testing"""
    
    profile = CRTMAtmosphericProfile{Float64}(
        n_levels = n_levels,
        pressure = zeros(n_levels),
        temperature = zeros(n_levels),
        water_vapor = zeros(n_levels),
        ozone = zeros(n_levels)
    )
    
    # Pressure levels (hPa, surface to TOA)
    p_surface = 1013.25
    for k in 1:n_levels
        profile.pressure[k] = p_surface * exp(-0.2 * (k-1)) * 100.0  # Convert to Pa
    end
    
    # Temperature profile (US Standard Atmosphere)
    for k in 1:n_levels
        p_hpa = profile.pressure[k] / 100.0
        if p_hpa > 226.32
            profile.temperature[k] = 288.15 - 6.5 * (p_surface - p_hpa) / 100.0
        else
            profile.temperature[k] = 216.65
        end
    end
    
    # Water vapor profile
    for k in 1:n_levels
        p_hpa = profile.pressure[k] / 100.0
        profile.water_vapor[k] = 0.015 * exp(-0.0001 * (p_surface - p_hpa))
    end
    
    # Ozone profile (simplified)
    for k in 1:n_levels
        p_hpa = profile.pressure[k] / 100.0
        if 10.0 <= p_hpa <= 100.0
            profile.ozone[k] = 8e-6 * exp(-((log(p_hpa) - log(30.0))/0.8)^2)
        else
            profile.ozone[k] = 1e-7
        end
    end
    
    return profile
end

# Surface type constants
const SURFACE_OCEAN = 1
const SURFACE_LAND = 2
const SURFACE_ICE = 3

function linear_index(i::Int, j::Int, k::Int, grid_size::Tuple{Int,Int,Int})
    nx, ny, nz = grid_size
    return (k-1) * nx * ny + (j-1) * nx + i
end

println("  ✓ All CRTM interface tests completed successfully")